Лабораторная работа 2

Сыромятников Дмитрий (КН-301)

Решить задачу Коши: $$ \begin{cases} \frac{dy}{dx} = 30y(x - 0.2)(x - 0.7) \\ y(0) = 0.1 \\ x \in [0, 1] \end{cases} $$

используя различные методы решения диффференциальных уравнений:

1. Явный метод Эйлера
2. Метод Эйлера с пересчетом
3. Метод Коши
4. Метод Рунге-Кутта 4-го порядка
5. Неявный метод Эйлера
6. Метод Тейлора 2 порядка точности на одном шаге
7. Метод Тейлора 3 порядка точности на одном шаге
8. 2-хшаговый метод Адамса
9. 3-хшаговый метод Адамса

Решения для каждого метода построить по 5, 10, 50 и 100 узлам.

Для применения методов Адамса осуществлять разгон методом Рунге-Кутта 4-го порядка.

Вычислить точное решение задачи. На графике отобразить ломанные, построенные перечисленными методами, точное решение. На основе графиков сделать выводы о погрешностях вычислений у перечисленных методов.

Решение

Вычисление точного решения

Общее решение

$\frac{dy}{dx} = 30 \cdot (x - 0.2)(x - 0.7)$ $| \cdot \frac{dx}{dy}$

$\frac{dy}{y} = 30(x - 0.2)(x - 0.7)dx$ $| \int$

$\int \frac{dy}{y} = \int 30(x - 0.2)(x - 0.7)dx$

$\ln |y| = 30(\frac{x^3}{3} - 0.9 \cdot \frac{x^2}{2} + 0.14 \cdot {x}) + C$

$\ln |y| = 10x^3 - 13.5x^2 + 4.2x + C$

$\Rightarrow y = C \cdot \exp (10x^3 - 13.5x^2 + 4.2x)$

Частное (искомое) решение

$y(0) = C \cdot 1 = 0.1 \rightarrow C = 0.1$

$\Rightarrow y(x) = 0.1 \cdot \exp (10x^3 - 13.5x^2 + 4.2x)$

from math import exp

def y_exact(x):
    """
    Точное решение дифференциального уравнения

    :param x: точка
    :return: значение решения в точке
    """
    return 0.1 * exp(10 * x ** 3 - 13.5 * x ** 2 + 4.2 * x)

Для вычисления решения методами Тейлора требуются частные производные правой части данного дифференциального уравнения ($f^{'}_x, f^{'}_y$):

$f(x, y) = 3y(x - 0.2)(x - 0.7)$

$f^{'}_x(x, y) = 6xy - 2.7y$

$f^{'}_y(x, y) = 3(x - 0.2)(x - 0.7)$

def f(x, y):
    """
    Функция f данного дифференциального уравнения
    """
    return 30 * y * (x - 0.2) * (x - 0.7)


def df_dx(x, y):
    """
    Частная производная функции f данного дифференциального уравнения по x
    """
    return 6 * x * y - 2.7 * y


def df_dy(x, y):
    """
    Частная производная функции f данного дифференциального уравнения по y
    """
    return 3 * (x - 0.2) * (x - 0.7)

Методы численного решения ОДУ

def _get_h(a, b, m):
    """
    :param a: начало отрезка
    :param b: конец отрезка
    :param m: число узлов по которым следует разбить отрезок
    :return: расстояние между соседними узлами
    """
    return (b - a) / (m - 1)


def euler_explicit(f, a, b, y0, n):
    """
    Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
    первого порядка явным методом Эйлера

    :param f: Функция правой части ОДУ вида dy/dx = f(x, y)
    :param a: Начальная точка интервала интегрирования
    :param b: Конечная точка интервала интегрирования (должна быть больше a)
    :param y0: Начальное условие (значение y в точке a)
    :param n: Количество узлов сетки на отрезке [a, b]
    :return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
    """
    h = _get_h(a, b, n)
    X = [a + i * h for i in range(n)]
    Y = [y0]
    for i in range(n - 1):
        Y.append(Y[i] + h * f(X[i], Y[i]))
    return X, Y


def euler_recalc(f, a, b, y0, n):
    """
    Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
    первого порядка методом Эйлера с пересчетом

    :param f: Функция правой части ОДУ вида dy/dx = f(x, y)
    :param a: Начальная точка интервала интегрирования
    :param b: Конечная точка интервала интегрирования (должна быть больше a)
    :param y0: Начальное условие (значение y в точке a)
    :param n: Количество узлов сетки на отрезке [a, b]
    :return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
    """
    X, Yt = euler_explicit(f, a, b, y0, n)
    h2 = _get_h(a, b, n) / 2
    Y = [y0]
    for i in range(1, n):
        Y.append(Y[i - 1] + h2 * (f(X[i - 1], Y[i - 1]) + f(X[i], Yt[i])))
    return X, Y


def cauchy(f, a, b, y0, n):
    """
    Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
    первого порядка явным методом Коши

    :param f: Функция правой части ОДУ вида dy/dx = f(x, y)
    :param a: Начальная точка интервала интегрирования
    :param b: Конечная точка интервала интегрирования (должна быть больше a)
    :param y0: Начальное условие (значение y в точке a)
    :param n: Количество узлов сетки на отрезке [a, b]
    :return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
    """
    X2, Y2 = euler_explicit(f, a, b, y0, 2 * n - 1)
    h = _get_h(a, b, n)
    Y = [y0]
    for i in range(1, n):
        k2 = 2 * i - 1
        Y.append(Y[i - 1] + h * f(X2[k2], Y2[k2]))
    return X2[::2], Y


def runge_kutta(f, a, b, y0, n):
    """
    Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
    первого порядка методом Ренге-Кутта 4-го порядка

    :param f: Функция правой части ОДУ вида dy/dx = f(x, y)
    :param a: Начальная точка интервала интегрирования
    :param b: Конечная точка интервала интегрирования (должна быть больше a)
    :param y0: Начальное условие (значение y в точке a)
    :param n: Количество узлов сетки на отрезке [a, b]
    :return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
    """
    h = _get_h(a, b, n)
    X = [a + h * i for i in range(n)]
    Y = [y0]
    for i in range(1, n):
        Xi = X[i - 1]
        Yi = Y[i - 1]
        K1 = f(Xi, Yi)
        K2 = f(Xi + h / 2, Yi + h * K1 / 2)
        K3 = f(Xi + h / 2, Yi + h * K2 / 2)
        K4 = f(Xi + h, Yi + h * K3)
        Y.append(Yi + (h / 6) * (K1 + 2 * K2 + 2 * K3 + K4))
    return X, Y


def euler_implicit(f, a, b, y0, n):
    """
    Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
    первого порядка неявным методом Эйлера

    :param f: Функция правой части ОДУ вида dy/dx = f(x, y)
    :param a: Начальная точка интервала интегрирования
    :param b: Конечная точка интервала интегрирования (должна быть больше a)
    :param y0: Начальное условие (значение y в точке a)
    :param n: Количество узлов сетки на отрезке [a, b]
    :return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
    """
    X, Ye = euler_explicit(f, a, b, y0, n)
    Y = [y0]
    h = _get_h(a, b, n)
    for i in range(1, n):
        Y.append(Y[i - 1] + h * f(X[i], Ye[i]))
    return X, Y


def taylor2(f, a, b, y0, n):
    """
    Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
    первого порядка методом Тейлора 2-го порядка

    :param f: Функция правой части ОДУ вида dy/dx = f(x, y)
    :param dfdx: Частная производная f по x
    :param dfdy: Частная производная f по y
    :param d2fdxdx: Вторая частная производная f по x
    :param d2fdydy: Вторая частная производная f по y
    :param d2fdxdy: Вторая смешанная производная f
    :param a: Начальная точка интервала интегрирования
    :param b: Конечная точка интервала интегрирования (должна быть больше a)
    :param y0: Начальное условие (значение y в точке a)
    :param n: Количество узлов сетки на отрезке [a, b]
    :return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
    """
    return euler_explicit(f, a, b, y0, n)


def taylor3(f, dfdx, dfdy, a, b, y0, n):
    """
    Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
    первого порядка методом Тейлора 3-го порядка

    :param f: Функция правой части ОДУ вида dy/dx = f(x, y)
    :param dfdx: частная производная f по x
    :param dfdy: частная производная f по y
    :param a: Начальная точка интервала интегрирования
    :param b: Конечная точка интервала интегрирования (должна быть больше a)
    :param y0: Начальное условие (значение y в точке a)
    :param n: Количество узлов сетки на отрезке [a, b]
    :return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
    """
    h = _get_h(a, b, n)
    X = [a + i * h for i in range(n)]
    Y = [y0]
    for i in range(0, n - 1):
        Xi = X[i]
        Yi = Y[i]
        y = Yi + h * f(Xi, Yi) + (h * h / 2) * (dfdx(Xi, Yi) + dfdy(Xi, Yi) * f(Xi, Yi))
        Y.append(y)
    return X, Y


def adams2(f, a, b, y0, n):
    """
    Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
    первого порядка методом 2-шаговым методом Адамса. Разгон производит, используя метод
    Рунге-Кутта 4-го порядка

    :param f: Функция правой части ОДУ вида dy/dx = f(x, y)
    :param a: Начальная точка интервала интегрирования
    :param b: Конечная точка интервала интегрирования (должна быть больше a)
    :param y0: Начальное условие (значение y в точке a)
    :param n: Количество узлов сетки на отрезке [a, b]
    :return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
    """
    h = _get_h(a, b, n)
    X = [a + i * h for i in range(n)]
    Y = [y0]
    X0 = X[0]
    Y0 = Y[0]
    K1 = f(X0, Y0)
    K2 = f(X0 + h / 2, Y0 + h * K1 / 2)
    K3 = f(X0 + h / 2, Y0 + h * K2 / 2)
    K4 = f(X0 + h, Y0 + h * K3)
    Y.append(Y0 + (h / 6) * (K1 + 2 * K2 + 2 * K3 + K4))

    for i in range(2, n):
        yi = Y[i - 1] + (h / 2) * (3 * f(X[i - 1], Y[i - 1]) - f(X[i - 2], Y[i - 2]))
        Y.append(yi)
    return X, Y


def adams3(f, a, b, y0, n):
    """
    Решает задачу Коши для обыкновенного дифференциального уравнения (ОДУ)
    первого порядка методом 3-шаговым методом Адамса. Разгон производит, используя метод
    Рунге-Кутта 4-го порядка

    :param f: Функция правой части ОДУ вида dy/dx = f(x, y)
    :param a: Начальная точка интервала интегрирования
    :param b: Конечная точка интервала интегрирования (должна быть больше a)
    :param y0: Начальное условие (значение y в точке a)
    :param n: Количество узлов сетки на отрезке [a, b]
    :return: Кортеж списков X, Y. X - узлы, Y - вычисленные значения функции в узлах
    """
    h = _get_h(a, b, n)
    X = [a + i * h for i in range(n)]
    Y = [y0]
    for i in range(2):
      Xi = X[i]
      Yi = Y[i]
      k1 = f(Xi, Yi)
      k2 = f(Xi + h/2, Yi + h * k1 / 2)
      k3 = f(Xi + h/2, Yi + h * k2 / 2)
      k4 = f(Xi + h, Yi + h * k3)
      y_next = Yi + (h / 6) * (k1 + 2*k2 + 2*k3 + k4)
      Y.append(y_next)

    for i in range(3, n):
        yi = Y[i - 1] + (h / 12) * (23 * f(X[i - 1], Y[i - 1]) - 16 * f(X[i - 2], Y[i - 2]) + 5 * f(X[i - 3], Y[i - 3]))
        Y.append(yi)
    return X, Y

Параметры

# Параметры задачи
a, b = 0, 1
y01 = 0.1
ns = [5, 10, 50, 100]
# ns = range(5, 100, 2) # важно начинать с 5 узлов
hs = [_get_h(a, b, n) for n in ns]

Данные

# Точные значения функции
X_exact = [0 + 0.001*i for i in range(1001)]
Y_exact = [y_exact(x) for x in X_exact]

# Расчет всех методов для всех n
data = {}
for n in ns:
    Xe, Ye = euler_explicit(f, a, b, y01, n)
    Xer, Yer = euler_recalc(f, a, b, y01, n)
    Xc, Yc = cauchy(f, a, b, y01, n)
    Xrc, Yrc = runge_kutta(f, a, b, y01, n)
    Xei, Yei = euler_implicit(f, a, b, y01, n)
    Xt2, Yt2 = taylor2(f, a, b, y01, n)
    Xt3, Yt3 = taylor3(f, df_dx, df_dy, a, b, y01, n)
    Xa2, Ya2 = adams2(f, a, b, y01, n)
    Xa3, Ya3 = adams3(f, a, b, y01, n)
    data[n] = {
        'euler_explicit': (Xe, Ye),
        'euler_recalc': (Xer, Yer),
        'cauchy': (Xc, Yc),
        'runge_kutta': (Xrc, Yrc),
        'euler_implicit': (Xei, Yei),
        'taylor2': (Xt2, Yt2),
        'taylor3': (Xt3, Yt3),
        'adams2': (Xa2, Ya2),
        'adams3': (Xa3, Ya3)
}

import plotly.graph_objects as go
from plotly.subplots import make_subplots

fig = make_subplots()

fig.add_trace(go.Scatter(
    x=X_exact, y=Y_exact,
    name='Точное решение',
    line=dict(width=4, color='#000000'),
    hovertemplate='%{y:.4f}<extra></extra>'
))

methods = [
    ('euler_explicit', 'Явный метод Эйлера', '#E6194B'),  # Ярко-красный
    ('euler_recalc', 'Метод Эйлера с пересчетом', '#3CB44B'),  # Зеленый "кричащий"
    ('cauchy', 'Метод Коши', '#FFE119'),  # Ярко-желтый
    ('runge_kutta', 'Метод Рунге-Кутта', '#4363D8'),  # Интенсивный синий
    ('euler_implicit', 'Неявный метод Эйлер', '#F58231'),  # Оранжевый
    ('taylor2', 'Метод Тейлора 2 порядка', '#00FF7F'),
    ('taylor3', 'Метод Тейлора 3 порядка', '#FFC0CB'),
    ('adams2', '2-шаговый метод Адамса', '#911EB4'),  # Фиолетовый
    ('adams3', '3-шаговый метод Адамса', '#42D4F4')  # Бирюзовый (циан)
]

for method_key, name, color in methods:
    X, Y = data[5][method_key]
    fig.add_trace(go.Scatter(
        x=X, y=Y,
        name=name,
        line=dict(width=2, color=color),
        hovertemplate='%{y:.4f}<extra></extra>'
    ))

frames = []
for n in ns:
    frame_data = []
    for method_key, name, color in methods:
        X, Y = data[n][method_key]
        frame_data.append(
            go.Scatter(x=X, y=Y, line=dict(color=color))
        )
    frames.append(go.Frame(
        data=frame_data,
        name=str(n),
        traces=[i for i in range(1, len(data[5]) + 1)]
    ))

fig.frames = frames

sliders = [{
    'active': 0,
    'currentvalue': {'prefix': 'Число узлов (n): '},
    'steps': [
        {
            'args': [
                [str(n)],
                {
                    'frame': {'duration': 0, 'redraw': True},
                    'mode': 'immediate',
                    'title.text': f'Сравнение методов (n={n}, h={h_val:.4f})'
                }
            ],
            'label': f'n={n} (h={h_val:.4f})',
            'method': 'animate'
        } for n, h_val in zip(ns, hs)
    ]
}]

# Обновление макета с динамическим заголовком
fig.update_layout(
    title=f'Сравнение методов решения ОДУ (n={ns[0]}, h={hs[0]:.4f})',
    xaxis_title='x',
    yaxis_title='y',
    hovermode='x unified',
    sliders=sliders,
    hoverlabel=dict(bgcolor='white', font_size=12)
)

fig.show()

Выводы

Построенные графики и их близость к точному решению полностью согласуются с погрешностями методов, с помощью которых эти графики решений построены: чем точнее метод, тем ближе решение, полученное с его применением к точному решению задачи Коши